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Abstract: The Karakoram Mountains are well known for their widespread surge-type glaciers and slight 
glacier mass gains. On the one hand, glaciers are one of the sensitive indicators of climate change, their 
area and thickness will adjust with climate change. On the other hand, glaciers provide freshwater 
resources for agricultural irrigation and hydroelectric generation in the downstream areas of the Shaksgam 
River Basin (SRB) in western China. The shrinkage of glaciers caused by climate change can significantly 
affect the security and sustainable development of regional water resources. In this study, we analyzed the 
changes in glacier area from 2000 to 2016 in the SRB using Landsat TM (Thematic Mapper)/ETM+ 
(Enhanced Mapper Plus) /OLI (Operational Land Imager) images. It is shown that the SRB contained 472 
glaciers, with an area of 1840.3 km?, in 2016. The glacier area decreased by 0.14%/a since 2000, and the 
shrinkage of glacier in the southeast, east and south directions were the most, while the northeast, north 
directions were the least. Debris-covered area accounted for 8.0% of the total glacier area. We estimated 
elevation and mass changes using the 1 arc-second SRTM (Shuttle Radar Topography Mission) DEM 
(Digital Elevation Model) (2000) and the resolution of 8 m HMA (High Mountain Asia) DEM (2016). An 
average thickness of 0.08 (+0.03) m/a, or a slight mass increase of 0.06 (+0.02) m we./a has been 
obtained since 2000. We found thinning was significantly lesser on the clean ice than the debris-covered 
ice. In addition, the elevation of glacier surface is spatially heterogeneous, showing that the accumulation 
of mass is dominant in high altitude regions, and the main mass loss is in low altitude regions, excluding 
the surge-type glacier. For surge-type glaciers, the mass may transfer from the reservoir to the receiving 
area rapidly when surges, then resulting in an advance of glacier terminus. The main surge mechanism is 
still unclear, it is worth noting that the surge did not increase the glacier mass in this study. 
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1 Introduction 


Glacier mass balance plays a prominent role in the study of the impact of climate change on 
glaciers and the response of glaciers to climate change. The local climate determines the supply 
and loss of mass and energy on the glacier surface and thus determines the mass balance of the 
glacier. Measurements of the glacier mass balance began in the Alps and Scandinavia in the 
1940s. In China, mass-balance measurements began in 1959 with the establishment of the 
Tianshan Glaciological Station, Chinese Academy of Sciences. To date, the continuous mass 
balance series, Urumqi Glacier No. 1, Qiyi Glacier and Dongkemadi Glacier, are longer than 20 
years in China (Wang et al., 2010; Yao et al., 2012; Wang et al., 2013, 2020). The traditional 
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glaciological method uses stakes or snow pits, which can obtain the mass balance at a single point 
by measuring the changes in the height of the stakes and the characteristics of the snow pits, and 
the point values can be extrapolated to the whole glacier to generate the glacier-wide mass 
balance (Pu et al., 2008). Glaciological mass balance is highly accurate; however, glaciers are 
usually located in uninhabited high mountains, and measuring stakes or snow pits is 
time-consuming and laborious, and the cost is too high. Therefore, the number of glaciers 
monitored by this method is very limited, and traditional methods have limitations in solving the 
glacier mass balance at the basin scale. With the gradual enrichment of remote sensing data for 
detecting glacial properties, research on glacier changes through remote sensing data has 
flourished. 

The Karakoram Mountains have a large ice volume, and glacier meltwater provides freshwater 
for many important inland rivers (e.g., the Yarkant River and Hetian River, both of which are the 
upper source of the Tarim River), which representing an important part of water resources in 
western China. According to the fifth report of the IPCC (Intergovernmental Panel on Climate 
Change), the average surface air temperature of the Earth has increased by 0.89°C over the past 
100 years (1901-2012) (IPCC, 2013). Climate warming accelerates the melting of glaciers, which 
not only threatens the production and life of local people but also affects the rise of sea level 
(Radić and Hock, 2011; Kääb et al., 2012). In the context of global warming, most of the world's 
glaciers have suffered mass losses. However, many studies have shown that the glaciers in 
Karakoram Mountains have experienced a balanced or slightly positive trend over the past decade 
(Gardelle et al., 2012a; Kääb et al., 2012; Berthier et al., 2019), which may be fostered stable or 
advancing glacier termini (Bolch et al., 2012; Rankl et al., 2014; Bhambri et al., 2017). This 
phenomenon is called the "Karakoram Anomaly" (Hewitt, 2005). 

Gardelle et al. (2012a) found that glacier mass balance of 0.11 (+0.22) m w.e./a for the 
Karakoram Mountains (glacier area is 5615.0 km’) during 2000-2008 using SRTM (Shuttle Radar 
Topography Mission) DEM (Digital Elevation Model) and SPOTS stereo image pair data 
(density: 900 kg/m? ). Kääb et al. (2012) further confirmed that glacier mass balance of —0.06 
(+0.04) m w.e./a (glacier area is 21,750.0 km’) during 2003-2008 by using ICESat laser altimetry 
(using a density of 900 kg/m? ), and Rankl and Braun (2016) proved that glacier mass balance of — 
0.08 (+0.10) m w.e./a (glacier area is 1107.0 km? ) for the period 2000-2012 using ice surface 
elevations of TanDEM-X and SRTM/X-SAR DEM (using a density of 850 kg/m? ). Another study 
using ICESat (Kaab et al., 2015) and ASTER DEMs from 2000 to 2016 (Brun et al., 2017) has 
shown that the Karakoram Mountains is located on the western edge of the West Kunlun 
Mountains, which is a larger mass-balance anomaly, the western Tibetan Plateau. These results 
indicate that glaciers in the Karakoram Mountains have experienced a stable state although the 
values slightly different each other. 

To overcome the scarcity of existing data sets, we updated the boundary data of glaciers in the 
SRB (Shaksgam River Basin) from 2000 to 2016. A geodetic method, based on a DEM difference 
of SRTM and HMA data, was used to determined changes in glacier mass balance during 
different periods. The function of debris cover in glacier mass changes was understood. Finally, 
combined with previous studies, discussed and analyzed the characteristics of surge-type glaciers 
and the climate factors of glacier changes. 


2 Study area 


The SRB (35°55'—36°49'N, 75°35'—77°30'E), on the northern slope of the Karakoram Mountains, 
is in western China (Fig. 1). The Shaksgam River, the main river in the SRB, originates from the 
Karakoram Mountains in Shenglinandaban and flows into the Yarkant River at Cha Hekou. The 
river course starts from the east-west direction, gradually changes to southeast after 30 km, and 
finally turns to the east-west direction in the last 40 km (Wang et al., 1989). This area is the upper 
source area of the Tarim River and is also known as "the Karakoram Corridor". The whole valley 
runs from southeast to northwest, with an average altitude of approximately 5100 m a.s.l. The 
world's second-highest peak, Mount Qogir (K2), is its main peak, with an altitude of 8611 m a.s.l. 
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The SRB is mainly controlled by a continental climate. According to the Chinese Second Glacier 
Inventory, the SRB contains 491 glaciers (in this paper, some small glaciers are removed) with a 
total area of 1884.2 km’. The most densely developed glacier area is Karakoram Mountains, 
China. Eight glaciers have areas less than 70.0 km’, which account for 55% of the total glacier 
area in the study area. The Yengisogat Glacier has an area of 359.1 km’ and a length of 42 km, 
and is the largest and longest glacier in China (Liu et al., 2015; Jiang et al., 2020). Due to the high 
altitude and low air temperature, measurement data are scarce. The Karakoram Mountains have a 
high altitude and steep terrain. Due to the interception of water and air by the high mountains and 
the role of storage in the mountains, there is more precipitation in the high-altitude areas than in 
the low altitude areas, which provides rich mass sources for the formation of glaciers and thus 
develops the largest and most concentrated glacier group in the middle latitudes (Xu et al., 2016). 
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Fig. 1 Location of the study area and distribution of glaciers in 2016. The background map is hillshade of 
SRTM (Shuttle Radar Topography Mission) DEM (Digital Elevation Model). The debris-covered ice was from 
Nico et al. (2018). 


3 Data 
3.1 Landsat data 


Landsat images (Table 1) were used to analyze glacier area changes in the SRB, including 
Landsat TM/ETM+/OLI. Eight Landsat TM/ETM+/OLI images with good quality were used in 
this study, which were downloaded from the website of USGS (United States Geological Survey; 
https://earthexplorer.usgs.gov). All scenes are from the USGS and are orthorectified with SRTM 
and ground control points from the GLS2005 (Global Land Survey 2005) data set, having 
achieved systematic radiometric and geometric accuracy. Due to the long-time span, high 
resolution and free download of Landsat data, it has been widely used in glacier remote sensing 
research (Paul, 2002; Bolch et al., 2010). 

3.2 Digital Elevation Models (DEMs) 


3.2.1 SRTM DEM 
The SRTM acquired interferometric synthetic aperture radar (InSAR) data simultaneously in both 


Table 1 Remote sensing images used in this study 


Satellite sensor Period of data Path/Row Cloud (%) Spatial resolution (m) 


Landsat TM 29 Aug 1998 148/35 4 30 
Landsat TM 4 Sep 2000 148/35 8 30 
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Landsat ETM+ 16 Jul 1999 148/35 13 15/30 
Landsat ETM+ 16 Jun 2000 148/35 2 15/30 
Landsat ETM+ 21 Jul 2001 148/35 6 15/30 
Landsat ETM+ 22 Jun 2002 148/35 3 15/30 
Landsat OLI 4 Jul 2015 148/35 3 15/30 
Landsat OLI 24 Sep 2016 148/35 6 15/30 


the C-band and X-band frequencies from 11 to 22 February 2000 (Farr et al., 2007). Data were 

acquired between 56°S and 60°N, accounting for 80% of the global land surface (Zyl, 2001). The 
1 arc-second (30 m) SRTM C-band, with linear vertical absolute height error of less than 16 m, 
and linear vertical relative height error of less than 10 m (at 90% confidence level). The SRTM 
DEM can be referred to the glacier surface in the last balance year (1999), with slight seasonal 
variation (Gardelle et al., 2013). The C-band exhibits some penetration of ice and snow, while the 
X-band exhibits less penetration. However, because the width of the X-band (50 km) is less than 
that of the C-band (225 km) (Farr et al., 2007), only 35% of the SRB glaciers are covered by 
X-band DEM. Low signal, signal-to-noise ratio etc. resulting in voids in the SRTM DEM data, 
but only 4% of the total glacier area was in the study. Hence, in this study, the 1 arc-second, 
non-void-filled SRTM C-band DEM was used to estimate the changes in glacier surface elevation 
after the 1999 ablation period. The data are available free from http://earthexplorer.usgs.gov/. 
3.2.2 HMA DEM 


The DEMs were generated from high-resolution along-and cross-track stereo images from 
DigitalGlobe satellites, with a resolution of 8 m, which were collected from 28 January 2002 to 
24 November 2016. The HMA DEM is composed of more than 4000 strips with a width of 13.1 
km and a length of 17.4 km, which reduces errors and eliminates seams. In this study, six images 
(100 kmx100 km) were collected during 2015-2016. HMA DEM can be referred to the glacier 
surface in the last balance year (2016), and the data were downloaded free from 
https://nsidc.org/datasHMA_DEM8m_MOS/Versions/1. Although HMA DEM suffers from data 
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Fig. 2 Resolution comparison of the two DEMs. (a), the 1 arc-second SRTM DEM from 2000; (b), the 8 m 


HMA (High Mountain Asia) DEM from 2015. The white spots in (b) are voids. The black line cycled areas were 
glaciers and both pictures were from Paul et al. (2019). 


4 Methodology 


4.1 Glacier boundary delineation 


Information on glacier boundaries is necessary to obtain changes in glacier mass balance. We 
selected the scenes during the ablation season (from June to September, mainly in August) to 
minimize the impact of snow cover and comprehensively considered the time, cloud proportion 
(Table 1). Seasonal snow and clouds can be distinguished by comparing images of the glacier 
body from neighboring years. 

Because of its simplicity and accuracy, the band ratio is the most suitable method to extract 
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glacier boundary (Albert, 2002; Paul et al., 2015). We use the method of band ratio threshold 
(TM3/TMS for TM/ETM+ and TM4/TM6 for OLI), with a threshold of 1.9 chosen for this task. 
The principle is to divide the surface object into glacierized and non-glacierized area by using the 
strong reflectivity of glaciers in the visible light band and the strong absorption in the 
near-infrared band. This method may mistake proglacial lakes for glaciers (Pan et al., 2012). 
Visual interpretation can easily correct this error. The surfaces of glaciers in the study area are 
covered by debris and supraglacial boulders, and the identification of these features is 
time-consuming and important work in later editing. There are several principles for 
distinguishing debris-covered glaciers. For example, Kääb (2005) suggested that areas with slopes 
lower than the threshold (23°) and directly connected to clean ice can be considered ice covered 
by debris. The colour of the surface moraine on the glacier is different from that of the debris on 
the outer edge of the glacier. The supraglacial boulders throughout a glacier tongue have the same 
spectral characteristics. With the increase in air temperature, the shape of the glacier tongue 
obviously decreases, which can be identified. Therefore, the above principles can be used to 
identify supraglacial boulders and debris-covered ice. At the same time, the process of extracting 
glacier information also refers to the second (Chinese glacier inventory 
(http://westdc.westgis.ac.cn/) and Google Earth. 

Although we have made extensive efforts to improve the accuracy of extracted glacier 
boundaries, errors still exist for various reasons. A certain degree of error assessment is needed. In 
this paper, the method proposed by Jin et al. (2005) is used to evaluate the error of the glacier 
boundaries. First, the buffer zone of the glacier boundary extracted from TM, ETM+ and OLI 
images are analyzed (taking half of the resolution of the remote sensing image as the buffer 
distance; thus, TM, ETM+ and OLI images have a buffer distance of 15 m). Then, the glacier area 
in the buffer zone is differentiated from the glacier area extracted from the original image, and the 
ratio between the result and the latter is the uncertainty of the glacier area extraction. The results 
showed that the uncertainty values of the extracted glacier boundaries in 2000 and 2016 were 
+3.72% and +3.54%, respectively, which were both within the range of previous research results 
(He and Yang, 2014). 


4.2 Geodetic mass balance calculations 


We use the geodetic method to obtain the glacier surface elevation changes by comparing the 
DEM data from different periods and then to estimate the glacier mass balance by combining the 
area and density of glaciers (Thibert et al., 2008; Fischer, 2011). 


M 
B=LY Ahs, (1) 
g i=l 
where B is the glacier mass balance (m w.e.); p is the ice density (kg/m’); S; is the glacier area 
(km®; M is the number of pixels contained in the glacier boundary; Ah; is the glacier surface 
elevation difference (m); and S, is the pixel size (m), that is, the spatial resolution of the DEM 
data. 


4.3 Co-registration and correction of DEMs 


Due to the differences in DEM acquisition, geodetic reference planes and processing methods, 
horizontal and vertical shifts might exist between the SRTM DEM and the HMA DEM. Any small 
horizontal shift will cause obvious DEM elevation errors in the mountain area. However, before 
estimating the elevation changes of the glacier surface, the relative registration of each DEM data 
is needed. Here, we used the method proposed by Nuth and Kääb (2011). Research shows that 
there is a cosine-curved obvious relationship between elevation difference, aspect and slope 
(Kaab, 2005) (Fig. 3a), based on the relationship relative vertical and horizontal distortions 
between the SRTM DEM and HMA DEM was corrected statistically (Nuth and Kaab, 2011). 
Based on the SRTM DEM, the spatial resolution of the HMA DEM was resampled to 30 m. 
Under the support of ArcGIS 10.2.2 software, a difference map was built with the SRTM C-band 
DEM and HMA DEM. On the basis of statistical analysis, 5% and 95% quantile thresholds are 
used to eliminate outliers DEM edges and near data gaps (Pieczonka and Bolch, 2015). Based on 
cosine fitting, we obtained the translation vector, and then the HMA DEM was moved to achieve 
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image registration. The calculation is as follows (Eqs. 2—6). In addition, based on the relationship 
between the elevation difference and the maximum curvature in non-glacier regions (Fig. 3b), the 
deviation, caused by the different spatial resolution of the two datasets, were refined for both 
glacier and non-glacier regions (Gardelle et al., 2012b). Before adjustment, the average elevation 
difference and standard deviation of the non-glacier regions were —19.83 and 14.91 m 
respectively. After these adjustments, the elevation difference in the non-glacier regions was 
focused on the average elevation difference at —0.38 m, and the standard deviation was reduced to 
11.8 m. After these corrections, the non-glacier regions had stabilized, which makes the DEMs 
suitable for estimating the change of glacier mass balance. 


dh 


oC aman ca, (2) 
fa OE, (3) 
tan (a) 


where dh is the individual elevation difference (m); a and ø represent the slope and aspect of the 
pixel (°), respectively; dh is the overall elevation bias between the SRTM DEM and HMA DEM 
(m), representing the vertical offset; a is the magnitude of the horizontal shift (m); b is the 
direction of the shift vector (°); c is the mean bias between the DEMs divided by the mean slope 
tangent (m), and ais the average slope of the DEM data (°). 

The offset in the x (m), y (m) and z (m) directions between the two DEMs can be calculated by 
the following formulas: 


x=axsin(b), (4) 
y=axcos(b), (5) 


z=cxtan(a), (6) 


where a, b and c can be obtained by cosine fitting. 
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Fig. 3 (a) Scatter plot of aspect vs slope for the standardized elevation differences for non-glacier areas. (b) 
Relationship between elevation difference and maximum curvature in the non-glacier areas of SRB (Shaksgam 
River Basin). 


4.4 Uncertainty assessment 


4.4.1 Evaluation of elevation change uncertainty 


The root mean square error or standard deviation of the elevation residuals between different 
DEMs in the area away from a glacier can estimate the uncertainty of the elevation changes 
(Bolch et al., 2011). 
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E = no glac (7) 


E= JE +E, (8) 


where E, is the average elevation bias of off glaciers (m); SDno glac is the standard deviation in the 
non-glacial areas (m); N is the number of measurement points; Em is the average elevation 
difference in non-glacier areas (m); and E is the error of elevation change (m). 

In order to calculate E, for a reasonable number of pixels near the glacier, we only used the 
pixels in the 5 km buffer zone around the glacier. Note that the DEM exhibits strong spatial 
autocorrelation in the elevation values of adjacent grids. To reduce the influence of 
autocorrelation, a decorrelation length based on the spatial resolution is suggested. For DEMs 
with a spatial resolution of 30 m, 600 m is suitable as the decorrelation length (Koblet et al., 
2010). 


Table 2 Statistical results of the vertical error before and after the adjustment of the HMA (High Mountain 
Asia) and SRTM (Shuttle Radar Topography Mission) DEMs (Digital Elevation Models) in non-glacier areas 


Before co-registration After co-registration 
N E, (m) E (m) 
En (m) SDno glac (m) En (m) SDao glac 
—19.83 14.91 —0.38 11.80 6960 0.14 0.41 


Note: Em, average elevation difference in non-glacier areas; SDyo giac, Standard deviation in the non-glacier areas; N, 
number of measurement points; E,, average elevation bias of off glaciers; E, error of elevation change. 


4.4.2 Effect of SRTM penetration 

The SRTM C-band (5.7 GHz) exhibits some penetration of ice and snow with a penetration depth 
between | and 10 m (Rignot et al., 2001). The penetration of the SRTM X-band (9.7 GHz) is 
expected to be low compared to that of the C-band (a hypothesis that still needs to be confirmed 
(Ulaby et al., 1986), we consider the elevation difference (SRTM C-band-—SRTM X-band) to be 
the first approximation of C-band radar penetration over glaciers (Gardelle et al., 2012b). Here, 
we estimate this penetration by differencing the SRTM C-band and X-band DEMs. The X-band 
and C-band are acquired simultaneously. There is no elevation change caused by ablation or 
accumulation (Gardelle et al., 2012b, 2013; Shangguan et al., 2015). Before estimating the 
penetration depth of the SRTM C-band, the SRTM X DEM and C DEM were co-registered and 
corrected. The final mean SRTM C-band penetration depth over the SRB was approximately 2.8 
m, which is consistent with the results of Gardelle et al. (2013), who reported a penetration depth 
of 3.4 m. 

4.4.3 Density assumption error 

The conversion from ice volume to mass is an important part when estimating glacier mass 
changes using the geodetic method (Huss, 2013). The glacier density must be considered to 
convert glacier elevation change to a mass balance change. Fresh snow, firn and ice exhibit large 
differences in density, and the density of glaciers varies, ranging from 100 to 917 kg/m? (Fischer, 
2011). The density parameter value will produce 5%—6% uncertainty in the glacier scale mass 
balance estimation results (Elsberg et al., 2001). In the Karakoram Mountains, it is assumed that 
the parameters used to convert the glacier elevation change into mass changes follow Sorge's law 
(Bader, 1954), that is, glacier density remains constant in the vertical direction, 900 kg/m” is used 
as the density conversion parameter, and the estimated mass balance was 0.11 (+0.22) m w.e. la 
during 1999-2008. If the glacial density in the accumulation zone is assumed to be 600 kg/m’, the 
mass balance drops to 0.05 (£0.16) m w.e./a, and both extreme density scenarios did not change 
the result of the proximity of the glacier area. In this study, a density of 850 kg/m’, with an 
uncertainty of 60 kg/m’, was used as the volume-to-mass conversion parameter. 


4.4.4 Accuracy of the glacier boundary 


The uncertainty of the glacier area refers to the error caused by the differences in S, and M in 
Equation 1, which is used to estimate the mass balance. To reduce this error, when calculating the 
changes in glacier volume, the maximum glacier area after the fusion of two periods can be taken; 
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when calculating the mass balance of glaciers, the union value of the glacier boundary from two 
periods can be taken (Gardelle et al., 2012b, 2013; Pieczonka et al., 2013). We finally obtained 
glacier areas in 2000 and 2016 of 1884.2 and 1840.3 km”, respectively. 

4.4.5 Uncertainty of mass balance 


The final mass balance uncertainty was determined using the estimated uncertainties of glacier 
area and elevation changes, and the glacier density uncertainty was 60 kg/m’: 


JE . 
j t Pe 


where Epps is the uncertainty of mass balance (m w.e./a); t is the observation period; p; is the ice 
density (850 kg/m’ ); Ap is the ice density uncertainty (60 kg/m? ); and p, is the water density 
(999.92 kg/m? ) (Pieczonka et al., 2013). 


5 Results 


5.1 Area changes 


There were 472 glaciers with a total area of 1840.3 km? in 2016 in the SRB (Fig. 4). In terms of 
area, large glaciers (=1 km? , occupy 94.3% of the total area) are dominant, while small glaciers 
(<1.0 km? , occupy 70.6% of the total number) are dominant in number (Fig. 4a). The glacier 
surface slope i in the SRB is approximately in the 16°—34° range, accounting for 95.6% of the total 
glacier area. Glaciers having a north, northeast, northwest or east aspects account for 63.2% of the 
total area. The shrank of the glacier on the southeast, east and south aspects were the most, while 
those on the west, northeast and north aspects were the least (Fig. 4b). 
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Fig. 4 Glacier distribution and change in the SRB. (a), area and number of 2016 glaciers in different size 
classes; (b), distribution of glacier area in different slope directions in 2000 and 2016; (c), the change of glacier 
area and area reduction rate under different size classes from 2000 to 2016; (d), the distribution of glacier area 
with elevation in 2016. 


The glaciers were significantly debris covered in the SRB and the distribution of 
debris-covered glaciers is shown in Figure 1. The total debris-covered area was 146.5 km’. 
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Among, all debris-covered glaciers, there were 20 glaciers with debris- covered areas larger than 
1.0 km’. Yengisogat Glacier has the most debris cover in the SRB (32.9 km’, about 9.3% of its 
area). About 86.2% of the debris-covered area is in the 4200 and 5200 m altitude range, 10.7% is 
higher than 5200 m, and only 3.1% is below 4200 m. 

Comparing the area of glaciers in 2000 with that in 2016, the reduction in glacier area was 43.9 
km’, with an annual average shrinkage rate of 0.14%/a. Small glaciers shrink significantly, 
however, the area loss of large glaciers is dominant in absolute area loss (Fig. 4c). Approximately 
91.5% of the total glacier area is found from 4400 to 6300 m, while only 1.6% is below 4400 m, 
and 6.2% above 6300 m. The median elevation is about 5682 m (Fig. 4d). The glacier topography 
analysis shows that the glacier area below 4000 m had disappeared completely, with an area of 
0.3 km’. Glaciers at low altitudes disappeared, and glaciers at high altitudes change very little and 
without disintegration, causing the decrease of the total number of glaciers. Glaciers in the SRB 
have experienced weakening area shrinkage from 2000 to 2016. 


5.2 Elevation and mass changes 


The glaciers in the SRB appear to show a weak surface thickening from 2000 to 2016. Glaciers, 
with an area of 1862.3 km’ (the union of the 2000 and 2016 glacier extents), experienced a 
thickening of 1.21 (40.41) m (0.08 (+0.03) m/a). Using 850 (+60) kg/m? as the density parameter, 
glacier mass gain of 0.06 (+0.02) m w.e./a. The distribution of glacier surface mass is uneven (Fig. 
5), which shows that the accumulation of mass is dominant in high altitude regions, and the loss of 
mass in low altitude regions (excluding the surge glacier). Yengisogat Glacier is typical. There are 
obvious glacier thickening phenomena in the upper reaches of the North ice flow (A in Fig. 5a), the 
upper reaches of the northwest ice flow (B in Fig. 5a), and the middle part of the main stream (C in 
Fig. 5a). 
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Fig. 5 Changes in glacier surface elevation in the SRB from 2000 to 2016. The glacier boundaries are based on 
the union of the 2000 and 2016 glacier area. These surge-type glacier and post-surge or quiescent during 1999— 
2011 were reported by Gardelle et al. (2012a, 2013). 


6 Discussion 


6.1 Comparison with previous study 


Previous studies agreed that the mass balance of glaciers in the Karakoram Mountains has 
maintained a weak increase or stable state since 2000 (Table. 3). Based on SRTM and SPOTS 
DEMs, a thickening of 0.12 (+0.19) m/a was found from 2000 to 2011 by Gardelle et al. (2013). 
Whereas a thinning of 0.09 (+0.12) m/a between 2000 and 2012 was reported over the Karakoram 
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by Rankl and Braun (2016), using Tan DEM_X and SRTM/X_ASR. And a thinning of 0.07 
(+0.04) m/a (2003-2008) was found in the Karakoram by Kääb et al. (2012), using SRTM and 
ICESat. In this study, SRTM and HMA DEMs acquisitions a mean thickening of 0.08 (+0.03) m/a 
between 2000 and 2016. This positive trend agrees with the results of Gardelle et al. (2013) but 
has slight differences from the results (a weak negative trend) of Kääb et al. (2012), and Rankl 
and Braun (2016). Different estimates of SRTM C-band penetration result in different thicknesses 
with those determined by Kaab et al. (2012). A mean SRTM C-band penetration of 2.4 (+0.3) m 
was used for the Karakoram Mountains (Kääb et al., 2012), less than 3.0 m (snow and ice), 8.0 m 
(accumulation areas) and an average penetration of 3.4 m in the central Karakoram Mountains 
(Gardelle et al., 2012a, 2013) and 2.8 m assumed in this study. The Karakoram Mountains is 
located at the west end of the Himalayan arc, stretching for thousands of kilometers, containing 
around 19,950.0 km? of glaciers (Cogley, 2011). The characteristics of glaciers in the eastern, 
western and central Karakoram Mountains are different. The study area of Gardelle et al. (2013) 
is mainly concentrated in the central Karakoram, while the study area of Kab et al. (2012) covers 
the whole Karakoram Mountains, our study area is mainly located at the eastern and central 
Karakoram Mountains. In addition, the differences in DEMs, the density of snow ice and glacier 
boundaries may explain part of these discrepancies. Previous studies estimated an average 
penetration depth of 2.7 m in the Karakoram Mountains (Zhou et al., 2017), 2.8 m in West 
Kunlun Mountains, 1.9 m in Pamir (Lin et al., 2017), 3.5 m in the Hindu Kush and 2.0-3.0 m in 
the Karakoram Mountains (Chen et al., 2018). In this study, the penetration depth was 2.8 m, in 
agreement with previous studies, which was considered acceptable and suitable for estimating 
glacier elevation changes in the SRB. However, these studies suggested that the average changes 
in glacier surface elevation of the Karakoram Mountains are certainly relatively small, even the 
glaciers have experienced a stable status in the Karakoram Mountains. 


Table 3 Comparison of mass balance budgets in this study and others 


Region Study period Elevation change (m/a) Mass balance (m w.e./a) Reference 
Karakoram Mountains 2003-2008 —0.07 (+0.04) —0.06 (+0.04) Kääb et al. (2012) 
Karakoram Mountains 2003-2009 —0.12 (+0.15) - Gardner et al. (2013) 

Central Karakoram 2000-2010 0.12 (+0.19) 0.10 (+0.16) Gardelle et al. (2013) 

Central Karakoram 2000-2012 —0.09 (+0.12) —0.08 (+0.10) Rankl and Braun (2016) 

Central Karakoram 1999-2008 - 0.11 (+0.22) Gardelle et al. (2012a) 
SRB 2000-2016 0.08 (+0.03) 0.06 (+0.02) This study 


Note: - indicates no data available. 


6.2 Debris-covered ice and clean ice 


The debris-covered areas account for about 8.0% of the total glacier (Fig. 6a). Ice cliffs, 
heterogeneous debris cover and supraglacial lakes could be the contributing factors to the 
complexity of the mass-loss patterns of debris-covered tongues (Pellicciotti et al., 2015). 
Although it is generally believed that a thick debris cover will act as a protective layer for glaciers 
due to its insulation (Benn et al., 2000), previous studies have reported that it will accelerate 
when the thickness of debris cover is smaller than a critical value (Ye et al., 2015; Zhang et al., 
2016). In this study, the debris cover has a positive influence on the ablation, and we found that 
melt was noticeably smaller on the clean ice than the debris-covered ice in the altitude 4000-5900 
m, during 2000-2016 (—0.63 (+0.03) vs 0.13 (+0.03) m/a) (Fig. 6b). Previous studies have found 
similar results in the Karakoram Mountains (Gardelle et al., 2012a), the eastern Pamir (Zhang et 
al., 2016), and the central Nyaingentanglha Range (Wu et al., 2019). 


6.3 Surge-type glaciers 


According to previous studies, at least 12 places surged in this basin (Fig. 5). Among the 12 
places, most surges occurred between 2000 and 2016 (Quincey et al., 2015; Bhambri et al., 2017), 
and one of them (C position in Figure 5a) occurred during 2011—2016. During 1999-2011, the C 
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position was defined with dashed black box in a post-surge or quiescent phase. Here, we found 
that an obvious thickening occurred at the C position, at the same time, there was obvious 
thinning in the upper reaches (below the A and B in Figure 5a), which is consistent with the 
characteristics of surge-type glaciers. In this study, we can't determine the specific time of the 
glacier surged here. However, according to the reported study, it is indicated that the C position 
(Fig. 5) of Yengisogat Glacier has a faster movement velocity of 900 m/a from June to August 
2012. Then, the glacier movement velocity decreased rapidly, that is, from June to August 2012, 
it may be at the peak of the surge phase (Quincey et al., 2015). 

Most of the glaciers either have stable termini locations or have been advancing, partly as a 
result of increased local precipitation (Janes et al., 2012) and partly due to glacier surging (Paul, 
2015) in the Karakoram Mountains. Here, the glacier surging did not increase the glacier mass. 
Generally, the rapid transfer of mass to the receiving area resulted in an advance of the glacier 
termini, but it was not always. Although some characteristics of the Karakoram glacier surges are 
similar to thermally controlled surges of other places (such as Svalbard), the main surge 
mechanism is still unclear (Quincey et al., 2015). 
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Fig. 6 Distribution of glacier area in 2016 (a) and glacier surface elevation changes during 2000-2016 (b) at 
100-m intervals by elevation in the SRB for debris-covered ice and clean ice 


6.4 Analysis of the factors of glacier anomaly 


Glaciers are a comprehensive product of climate, topography and other factors. Since the terrain 
has not changed in recent decades or even hundreds of years, glaciers are sensitive to climate 
change (Shi, 2000). Some potential explanations, such as cloud cover, precipitation and 
temperature can be found in the current studies, which may be driving factors of the Karakoram 
glacier changes. Based on the global meteorological reanalysis data and meteorological station 
data, Forsythe et al. (2015) and Bashir et al. (2017) have reported that an increase of cloudiness in 
the Karakoram Mountains. The increase of cloud cover may reduce the incoming shortwave 
radiation (Forsythe et al., 2015). However, incoming shortwave radiation is generally considered 
to be the main energy source of glacier mass changes (Forsythe et al., 2015). In addition, previous 
studies have reported that increased cloud cover leads to increase in precipitation (Forsythe et al., 
2015; Bashir et al., 2017). The climate of Karakoram is mainly influenced by westerly circulation 
and Tibetan anticyclone (Archer and Fowler, 2004). Cannon et al. (2015) reported that the 
frequency and intensity of such precipitation by westerlies-dominated have increased in recent 
years, which seems to lead to a slight increase in winter snowfall in Karakoram (Norris et al., 
2019). In turn, precipitation change is considered to be a powerful control on glacier mass balance 
in the Karakoram Mountains (Kapnick et al., 2014; Forsythe et al., 2017). Moreover, in recent 
decades, summer cooling has been a particularly important driving factor for the mass balance of 
glaciers in the Karakoram Mountains (Fowler et al., 2010; Forsythe et al., 2017). It is worth 
noting that the tree-ring chronology did not provide any evidence that the temperature of the 
Karakoram Mountains is not consistent with the rest regions in High Mountain Asia. 
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7 Conclusions 


This paper takes the glaciers in the SRB as the research object, and the glacier boundaries in 2000 
and 2016 were extracted according to Landsat TM/ETM+/OLI images. Based on HMA DEM and 
SRTM DEM data, changes in surface elevation and mass balance in this study have been 
estimated from 2000 to 2016. 

We found that the SRB contained 472 glaciers, with a total glacier area of 1840.3 km? in 2016. 
The reduction in glacier area was 43.9 km? with an annual average shrinkage rate of 0.14%/a, and 
glacier area mainly decreased in the southeast, east and south directions, during 2000-2016. 
Compared with the melt of mountain glaciers in western China, even global, the glaciers in the 
SRB have experienced a slight shrinkage. In addition, the debris-covered area was 146.5 km’, 
accounting for about 8.0% of the total area. 

The glacier thickness showed a weak increasing trend, with a thickening value of 0.08 (+0.03) 
m/a, or 0.06 (+0.02) m w.e./a. The thinning value was significantly greater on the debris-covered 
ice than the clean ice in the altitude 4000-5900 m, between 2000 and 2016 in the SRB (—0.63 
(+0.03) vs 0.13 (+0.03) m/a). Research result shows that the glacier mass was rapidly transferring 
from the reservoir to the receiving area, resulting in an advance of the glacier termini, but not 
always. 
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